Multistable behavior above synchronization in a locally coupled Kuramoto model 
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A system of nearest neighbors Kuramoto-like coupled oscillators placed in a ring is studied above 
the critical synchronization transition. We find a richness of solutions when the coupling increases, 
which exists only within a solvability region (SR). We also find that they posses different character- 
istics, depending on the section of the boundary of the SR where the solutions appear. We study 
the birth of these solutions and how they evolve when K increases, and determine the diagram of 
solutions in phase space. 
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I. INTRODUCTION 



The ubiquity of phenomena linked to coupled chaotic 
systems has made them the focus of interest for the last 
twenty years. The study of these systems has raised inter- 
est with the intent of realistically modeling spatially ex- 
tended systems, as diverse as Josephson junction arrays, 
multimode lasers, vortex dynamics, biological informa- 
tion processing, neurodynamics as well as applications 
in communications ^{5], with the belief that dominant 
features will be retained in such simple models. Coupled 
systems with local interactions are of special importance. 
In particular the Kuramoto model[n] in its local version: 
locally coupled Kuramoto model (LCKM) has raised at- 
tention since most features of systems with phase cou- 
pling appear in this particularly simple model [TDHIl]- 

This system which shows synchronization in the mean 
frequency has been thoroughly studied at and before full 
synchronization. The oscillators cluster in average fre- 
quency decreasing the number of clusters until they come 
into a single cluster at full synchronization. At this mo- 
ment the frequency becomes constant and the phases 
lock, such that all phase differences are constant. The so- 
lution for full synchronization and its stability has been 
studied by many authors [TSHTO] . Zheng et al [lT7 al- 
ready in 1998 pointed out that the behavior of the order 
parameter " indicates the coexistence of multiple attrac- 
tors of phase locking states" above the synchronization 
critical coupling, Ks- Synchronization with coexistence 
of attractors has been reported by different authors and 
in different fields below ptiT4l EUS] but, to the 
best of our knowledge nobody else has pursued the mat- 
ter above full synchronization. If we have the intention 
to simulate real systems, mostly in technological applica- 
tions, it is important to know whether or not we can move 
freely below and also above synchronization within stable 
solutions and to know whether or not they are unique. 
In this work we shall study a LCKM above complete syn- 
chronization, we shall show that there is an unexpected 
richness of behavior: multistable solutions appear and 



we cannot change the strength of the coupling without 
danger of falling into different attractors. 



II. LOCALLY COUPLED KURAMOTO MODEL 
ON THE SYNCHRONIZED REGION 

The model that we use is described by the following 
equations: 



K [sin ( 
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where n = 1, A^, and w„ G {cj} is the set of natural 
frequencies of the oscillators. The ring topology is defined 
by the periodic conditions ^at+i = 9i and — (^n ■ There 
is a minimum value for the coupling constant K, denoted 
as critical synchronization coupling Ks, that drives the 
system into a fully synchronized state 16 18 . In this 
state the oscillators instantaneous frequencies assume a 
constant value D, = SjLi that remains unchanged 
for any K e [Ks,oo). 

The set of equations ([!]) in a synchronized state may 
be written as 
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with <j)n ~ On — On+1 and n — 1,...,N — 1 that en- 
ables one to write the phase locking conditions as (/)„ — 
(j)n {K , {uj}) , with multiple stable solutions depending 
both on the number of oscillators A^ and on the coupling 
constant K > Kg. On a synchronized state any 0„ may 
be written as a function of an arbitrarily chosen distance 
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for n — 1, ...,n* — 1 and 
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for n = n* + 1, N-l. With the identity J^f^^ = 
it is possible to write 4>m ss the sum of all others (pm and 
the structure of the system allows us to reduce the set 
of A'' equations ([2| to a single equation on two variables 
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The choice of 0„. will become apparent in the next para- 
graph. 

If one removes a single link (interaction) between any 
pair of oscillators the result is a chain with free ends. 
If we can find the chain encapsulated inside the ring 
with the highest synchronization coupling we may be 
able to find the solutions above Kg. For each of the 
N possible ways for which this procedure may be done, 
we label the oscillators as wi, wat, keeping their or- 
der when we change the labels. Following the calculation 
of Strogatz and Mirollo [TS] we find the coupling con- 
stant at complete synchronization for a given chain as 



Ks{l) = maxi<i<Ar 



Therefore we iden- 



tify n* as the value of I when this expression is the high- 
est among all the synchronization couplings Ks{l) and 
we call it i^f 

We write equation Q as follows: 
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where s 



and solve it numeri- 



cally using Mathematica. The general form of the solu- 
tion is presented in figure [l]for iV = 10, 20 and 50 where 
the frequencies {wj^y were generated from a uniform dis- 
tribution function defined in the interval [—10,10]. We 
have performed calculations for different natural frequen- 
cies realizations, all the results present the same features 
as will be described, therefore we keep here to a single 
realization for a given TV to avoid confusion. 

We see that there are multiple phase locking solutions 
for the system above Kg being spontaneously generated 
on a confined region on the bifurcation diagram (BD), 
as indicated by the dotted lines. These solutions are of 
two kinds: a) we call type I solutions the ones generated 
on the bottom boundary of BD with sin ^ s that 
bifurcate into branches keeping sign (cos <j)n*) invariant; 
b) solutions with opposite signs of cos 4>n* sharing the 
same origin on the top of BD (close to sin = s) are 
called type II. 

The confinement region may be obtained realizing that 
for each value of sin (/>„» there is a maximum value of K 
below which ([s]) is never satisfied. By assuming that 
sin(/)7v = — s we find that all solutions will appear at 
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The equal sign describes the points where type I solu- 
tions touch this bordering line and is represented by the 
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FIG. 1: (Color online) Numerical solutions of Q representing 
the bifurcation diagram for a ring of oscillators (stability not 
explicit): Top: A*' = 10 with natural frequencies cji = 6.9, 
u)2 ~ 2.8, ui-i = —0.4, CJ4 = —2.6, CJ5 = 1.3, ujq — —6.8, 
UJ7 — 0.8, ojs = —1.6, iug = —9.5 and cjio = —6.7. Middle: 
realization with A'^ = 20. Bottom: realization with A'^ = 50. 
The sets of natural frequencies {(^jjv are generated from a 
uniform distribution function defined in the interval [—10, 10]. 
Black lines represent solutions with cos (j)n* > 0, blue lines 
are have cos <f)n* < and the dotted lines are the limiting 
boundaries sincijv = —s and sin(A„* = s. 



inferior dotted line in BD. The solvability region (SR) 
contains all the synchronized solutions and it can be de- 
fined by Q. It is worth mentioning that equation ^ 
shows that the critical synchronization coupling for the 

^chain 

ring satisfies the condition Ks > — . 

A closer look at the BD will show that there is always 
one solution from each bifurcation that is tangent the 
boundary of the SR (which we call SB) at only one point 
(this will appear clearly when we discuss figure [5|: it cor- 
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responds to sin (/17V = —s for type I and sin(/)„* = s for 
type II solutions. On all of these tangent points, inde- 
pendent of the boundary, a straightforward calculation 
will show that they satisfy the condition 



cos 
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The number of bifurcations depends on the number of so- 
lutions of equation ([t]) over each boundary of SR. Since 
the cosine argument may be expanded as a Laurent series 

defined by (l)n — ^0 + J2m=i decaying behavior 

as a function of K guarantees a finite number of solu- 
tions. From the definition of the phase differences ([s]) it 
is possible to see that the ^0 term depends on the size 
of the system so that the effect of increasing the number 
of oscillators also increases the number of phase locking 
solutions above Kg, as may be observed in figure [l] 

One way of addressing the question of how the multiple 
solutions are generated is to look at the solutions of 0„* 
and 0JV on the tangent points. Starting with the type 
I solutions the condition sini/ijv — —s imposes that 
should satisfy the equation 



4^)(^,) = -.s ^'"^^ + ^^\ m,.0,l,2,... (8) 

For each possible value of mi equation ([S]) admits two 
solutions: 
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If we fix (pn* = 0!/.'^'' then the values assumed by mi 
will provide all the values of K where the solutions with 
cos 0„* > are tangent to the limiting boundary. Simi- 
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larly when 



the values of mi give the tangent 



points for the solutions with cos0„. < 0. Figure [2] (top) 
shows these solutions for the {i^}jv=io realization used 
previously, when = 10. 

The same procedure may be performed for the type II 
solutions, where the condition sin — s imposes 

(i)\.'{m2)^-S ^ , 7712=0,1,2,... (10) 

Now for each value of TO2 we also have two types of so- 
lutions for (piq: 
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The values of K for the tangent points on this bound- 
ary are also shown on figure [2] (bottom), for the same 
realization with N = 10. 
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FIG. 2: (Color online) Graphical representation of the nu- 
merical solutions of I?! for the A'' = 10 realization, with 

y = s 1) (blue curve) and fiji = cos {Y^'^I^^, (pn^ 

(black curve) calculated on the tangent points from type I and 



II solutions, respectively. Top: 



solutions (black dots) at 



K(^^^o) = 9.9905, A'(„i=i) = 23.0467, A'(„i=2) = 190.016 



and 



solutions (blue dots) at 



(mi=0) 



13.443 and 



(ii.-i-) 

N 

10.75473, -?s:{„i2=i) = 84.3981 and 



iS'(„j^i) = 104.7171. Bottom 



solutions (black dots) 



at J('(m2=o) 
tions (blue dot) at K^^^^2) = 31.1029 



solu- 



With the description of all multiple phase locking so- 
lutions above (at least on the tangent points) it is 
possible to have some insight of their origin. If we go 
back to equation Q and look at the first term on the left 
hand side it is possible to visualize that as K increases 
the argument goes beyond 27r but not simply adding to 
it a multiple of 27r. In this way it is the presence of this 
term - connecting the first to the last oscillators of the 
chain - that generates the multiple solutions. The no- 
table symmetry between the type I and type II tangent 
point solutions comes from the fact that if we remove the 
interaction term between oscillators 9i and On the criti- 
cal synchronization coupling of the resulting chain is the 
same as the one obtained from extracting the interaction 
term connecting oscillators On* and 0„*+i, even when the 
phase locked solutions may be different. Nevertheless this 
symmetry is not perfect because the number of type I and 
type II solutions are in general not the same. 
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III. STABILITY OF SOLUTIONS AND BASINS 
OF ATTRACTION 



To perform a linear stability analysis of the solutions 
it is necessary to obtain the jacobian matrix, but as the 
equations of motion are invariant by global phase trans- 
lations 9,1 — > ^'ji + (for every n) the analysis is a little 
more complicated. 

It is necessary to realize that the freedom of gauge 
reduces the N equations of motion ([ij to a — 1 di- 
mensional system (this is an effect of the constraint 
X]n=i ^1 = 0)- Although a gauge fixing condition would 
enable us to eliminate the extra degree of freedom and 
perform the stability analysis on the On variables we be- 
lieve there is a more tractable way to do this. Instead of 
considering the equations of motion as presented in ([T]) 
we write the equations on the 0„ variables and eliminate 
the degree of freedom replacing ^at by —J2n=i 't'ri- As 
only the equations for i/ii and (pN-i are dependent on 
(^AT, we have 



and for n = 1 . . . . , — 3 



—K sin I ^ (j)n I + 2 sin 01 — sin 02 , (12a) 

\n=l / 

-\-K sin j (pn I + sin (/)jv-2 ^ 2 sin (/)jv-i ■(12b) 

\n=l / 

The rest of the equations for n = 2,...,Af — 2 do not 
depend on (/)jv and are expressed in the usual form 

+ K (sin (j)n-i — 2 sin (j)n -f sin (j)n+i) ■ (12c) 

Due to the structure of the equations, the elements of 
the jacobian matrix may be written as Jn,m — d(f,^(l)n- 
Since most of the elements are null let us focus on the 
nonzero ones: the first and last lines are complete with 
the elements given by 
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Jl,2 — K cos 01 — COS ( ^ 



Jn-i.n-2 = K 




Jn-1, 



N-l 
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2 COS 0Ar_i -I- COS ^ (f^j 



(13a) 
(13b) 
(13c) 
,(13d) 



Jl.n 
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all other lines from n = 2, ...,N—2 have nonzero elements 
only on the diagonal and first neighbors, namely 

Jn,m — K COS 4>rn {^m,n-l — '^^m,n + ^m,n+l) ■ (13f) 

The linear stability of the fixed points is determined by 
the structure of the jacobian eigenvalues: a solution will 
be stable only if the real part of all eigenvalues are neg- 
ative. Since a general system with N oscillators has an 
{N — 1) X (N — 1) jacobian matrix, the eigenvalue equa- 
tion is too complicated to be treated analytically. In 
this way we performed a numerical approach that con- 
sists in fixing a value for the coupling K and replacing 
the phase locked solutions 0„ on the jacobian matrix ele- 
ments by the numerical calculation of the eigenvalues A„, 
for n = 1, iV — 1. Figure |3] shows the real part of the 
largest eigenvalue for all fixed points for the case = 10 
with K = 120 (the fixed points may be inferred either 
from figure [l]or fi gure El) . Repeating the procedure for 
other values of K on the range [Kg, 120] we were able to 
determine the stability of all solutions, as shown on the 
BD from figure |4] Our results confirm that the stability 
of a branch is not altered by variations of K, therefore 
the results obtained from figure |3] may be extrapolated 
to the whole BD. 
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FIG. 3: (Color online) Largest jacobian eigenvalue for the 
A'' = 10 case with K = 120. Black dots represent phase lock- 
ing solutions with cos (j)n* > while blue squares represent 
phase locking solutions with cos<^„» < 0. 

Let us focus on the A^ = 10 case. If we look at the 
type I solutions (boundary of the SR with sin0„. ^ s), 
as shown on top of figure|4] we observe that if cos 0„* > 0, 
two branches are born on a saddle node bifurcation and 
the stable branch is tangent to the SB curve, sin 0jv = 
(details are shown in Fig. 5 left). On the other hand both 
branches with cos 0„* < born at the same point, are 
unstable with at least one positive eigenvalue. If we take 
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FIG. 4: (Color online) Bifurcation diagram for A'' = 10 case 
with explicit stability. Continuous (dashed) curves represent 
stable (unstable) solutions. Top; type I solutions. Bottom: 
type II solutions. 



random distributed natural frequencies on a ring (as 
shown in [17)). which makes it impossible to have a bi- 
furcation with a sine equal to ±1 Nevertheless there 
is a sine equal to ±1 near the bifurcation, with the dif- 
ference becoming smaller as — >■ cx3 in agreement with 
previous literature where this fact has played a crucial 
role [ini US] . On the next section we will show how these 
deviations and the multiple solutions are generated from 
a chain of oscillators. 
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FIG. 5: (Color online) Zoom of the first two bifurcations 
present on the BD of the A'^ = 10 case. Left: Zoom of the 
Ks bifurcation from figure |4] (top) with a tangent point from 
the stable solution touching the SR boundary at sin <j)M = —s. 
Right: first bifurcation from type II solutions from figure |4] 
(bottom) showing that the stable solution is actually born 
with cos (j)n* < and changes sign when it touches the SR 
curve sin (An* = s. 




a look at the bottom of figure [4] it is possible to see that 
the type II solutions (close to sin(/)„* — s) also present 
bifurcations that generate two unstable solutions. The 
stable solutions also have a point tangent to a SR curve 
(sin0„. = s), but since sin0„» = s, cos(/)„* changes sign 
at the tangent point, different from the type I case, the 
stable solution starts at the bifurcation with cos 0„. < 
and changes sign at the sin — s line. 

Now lets analyze the solutions near the bifurcations. 
All branches are born at local minima of the function 
K ((pn') defined in the synchronized region (0„ = Q), 
where the critical synchronization coupling is the ab- 
solute minimum. Bearing this in mind, when we take a 
closer look at bifurcations from each type of solution an 
interesting feature is observed (figure [s]) : sin 0„ is not 
equal to ±1 at Kg or at any other local minima for any 
n. The explanation for this behavior lies on equation 
([5]): in order to have a bifurcation with a sine equal to 
±1 it would necessarily be either 0„. or (j)]s[ to present 
this property, but since 0Ar is a nonlinear function of 
there are accessible synchronized solutions prior to the 
appearance of the sine equal to ±1. This result is in 
agreement with the condition 
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The general picture we obtained for the BD with ran- 
dom distributed frequencies, from both simulation and 
numerical calculation, may be summarized as follows: for 
a given N, a configuration {w}^ determines a SR for the 
system in which all solutions come from bifurcations near 
the solvability boundaries; all stable solutions correspond 
to branches that touch the solvability boundary; it is pos- 
sible to extrapolate the definition of the two types of so- 
lutions to the bifurcations themselves and label then as 
Kj, where I — I 01 II, for the type I or type II bifurca- 
tions and j = 1, 2, gives the order of the minima; odd 
values of j correspond to saddle node bifurcations while 
even numbered bifurcations generate two unstable solu- 
tions; none of the bifurcations may appear with a sine 



equal to one because condition ( 14 ) must always be sat 



found to be satisfied at the critical coupling Ks by any 



isfied; in addition to the stability analysis our numerical 
simulations showed that every stable tangent point sat- 
isfy the condition Ok cos (^n^n- '^n) > although we 
are not able to explain why it happens at this moment. 

With the stability of the solutions fully described we 
turn our attention to the basin of attraction of the sta- 
ble solutions. Since a general LCKM is a high dimen- 
sional system a usual approach consisting of a graphical 
analysis becomes extremely difficult either for numerical 
computation or graphical visualization. To outline these 
difficulties we consider a statistical approach: given that 
the fixed points are represented by values of sin (/)„» we 
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generate a large sample of random initial conditions for 
the phases on the interval [— 7r,7r], for fixed values of K, 
and estimate the size of the basin of attraction by the 
probability of the system to reach each of the stable so- 
lutions. 

Starting with the iV = 10 case we could observe that as 
long as the values assumed by K lie on a region where the 
number of stable solutions is kept constant, the size of the 
basins of attraction presents only small statistical fluctu- 
ations, which lead us to believe that these variations of K 
have little effect on the size of a basin. When we vary the 
coupling constant across regions with increasing number 
of solutions we observe that when a stable fixed point is 
created its basin of attraction steals the majority of its 
size from the closest fixed point (in the sin(/)„* space), as 
it is illustrated on figure [6] for some values of K. 



Pr(sin 4>n') 




FIG. 7: (Color online) Probability of the system reaching each 
of the stable fixed points (for large values of K) representing 
the estimation of the relative basin size. Blue dots represent 
A'' = 50, lilac squares represent A'^ — 100 and beige diamonds 
represent A'^ = 200. 
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FIG. 6: (Color online) Estimation of the basin size for each 
of the stable fixed points within the A'' = 10 case. Blue dots 
represent K = 20, lilac squares represent K — 50 and beige 
diamonds represent K = 120. 



We found numerically that for a given system the size 
of the basins of attraction is not evenly distributed among 
all fixed points: the stable solutions with phase locking 
sin (pn* going to zero as JsT — )■ oo attract the majority of 
the initial conditions. Fortunately the behavior of the 
system for large values of K is easier to analyze: in the 
limit _ftr — >■ cx) we have 0„ = (f>, for n — I, A'^ — 1, with 
(j) given by the solutions of 



sin [{N - 1) (/)] 



(15) 



coming in two types. 
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In this regime the solutions depend only on the system 
size A^ and the values assumed by the probability of 
reaching a stable fixed point are within the curve defined 
by / (sin0„*) ~ exp (— asin^ as shown on figure [t] 



IV. CONNECTION TO THE CHAIN 

In the previous sections we described the synchronized 
region of a ring of oscillators and showed how some of the 
properties of the chain are still present in this topology. 
Now we focus our attention to understand the appear- 
ance of multiple solutions starting from the solution of a 
chain. 

To construct the chain we remove the link between os- 
cillators 1 and N in the ring to end up with a chain of 
oscillators where the critical synchronization is defined 

by sin0„-» = s^j^ — (note that the sin^^r term does not 
enter in the dynamics). A saddle node bifurcation ap- 
pears at K — X^hain g^^^ naturally we have sin0„. = s. 
The effect of closing the chain into a ring by connecting 
these two oscillators inserts the extra term present on 
the left hand side of (5]), which is a nonlinear function of 
. Instead of a single solution for the fully synchronized 
state that extends for all K > Kg, this new configura- 
tion generates multiple stable solutions born at the local 
minima of the implicit multivalued function K {(pn* ) that 
bifurcate into pairs of solutions as K increases. 

Now we build the ring from the open chain in a con- 
trolled way by coupling a continuous parameter a e [0, 1] 
to the interaction term sin m equation ([2| . It is easy 

to see for a chain (a = 0) that sin0„. = — for 
all K > Kf^'^™ [TS]. The effect of increasing a deforms 
the original solution (as shown in figure |8]) which in turn 
will generate the others observed for the ring. For small 
values of a a hysteresis figure appears at a fold bifur- 
cation (fig. [s^), creating the K^^'' bifurcation from the 

stable branch born at K^^^^ , but as a increases the turn- 
ing point goes to infinity (fig. [8]d), completely separating 

the original into two stable solutions. Bifurcations ^ 

and ifg^' do not present such turning points, as they seem 
to come from K ^ oo (figs, [sja andjs]:). 

Higher values of a give birth to a small closed circuit 
near = 7r/2 ([s]; and [9]), which via a deformation 
process increases the circuit until the turning points go 
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FIG. 8: (Color online) Bifurcation diagrams ^n* versus K 
for the A'^ = 10 case as a function of the enclosing parameter 
a. (a) Hysteresis figure showing the cusp catastrophe when 
the stable solution loses stability and disappears; K2 coming 
from infinity {a = 0.3). (b) Turning point goes to infinity 
(bottom) creating two distinct solutions and comes from 
infinity (a — 0.5). (c) Birth of A'l^ and K^' from the closed 
circuit (a = 0.7). (d) Bifurcation diagram of the ring topology 
showing all solutions on the region (a = 1.0). 
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FIG. 9: (Color online) Closer look at the region of the birth 
of A"!^ and K^^ showing solutions for four different values of 
a: ai = 0.683, Q2 = 0.7, 03 = 0.75 and Q4 = 0.8. 

to K — >■ 00, creating K2^^ and ^^3^^^ each having two 
solutions. Figure [9] shows how the closed circuit of un- 
stable solutions is born {a — ai) and as the parameter 
increases the solution crosses the 0„. = 7r/2 line, creat- 
ing one stable solution; at this point the bubble starts to 
deform creating the turning points on both sides of 7r/2 
and generating four solutions. When a — 1 the turn- 
ing points go to infinity and the system presents all the 
properties of the ring previously described. 

The complete structure of the synchronized state as 
the system is driven from a chain into a ring may be 
visualized on the stability diagram shown in figure |10[ 
small values of a change the phase locking solutions (/)„ 
but do not alter the main structure of the chain; for a ~ 
0.1 and iiT ~ 35 (region A) the properties of the ring start 
to become apparent as the hysteresis appears and the 
system presents bi-stability on a closed interval of values 
of K; regions Bl and B2 show the bifurcations and 

K^^^ coming from infinity and the disappearance of the 
hysteresis; region C shows the birth of the closed circuit 
that creates if 2^^' and K^^^\ 

V. CONCLUSIONS 

In conclusion we studied a locally coupled Kuramoto 
model above the full synchronization transition for a ring 
of chaotic oscillators. We found a very rich panorama of 
solutions, although there is only one at the critical value 
for full synchronization. We were able to determine (an- 
alytically) the solvability region (SR) where the solutions 
exist and to show (numerically) that they all come from 
minima of the function K (0„. ) , with being a spe- 
cific chosen phase difference. From the observation that 
every bifurcation has a solution branch with a point tan- 
gent to the SR curves we were able to show that the 
multiplicity of solutions comes as distinct discrete values 
assumed by and ipN- 

The stability analysis of the solutions showed the exis- 
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FIG. 10: (Color online) Stability diagram for the system. The 
birth of a pair of solutions is represented by a continuous 
line as one goes from lower to higher values of K, while the 
destruction of a pair of solutions is represented by dotted 
lines (the definitions are inverted if one goes from right to 
left). Thick (black) lines denote the creation of one stable 
solution and the thin (blue) ones account for the bifurcations 
with no stable equilibria. A: birth of the hysteresis. Bl and 
B2: bifurcations that come from infinity. C: Birth of the 
closed circuit that generates K2^^^ and k!^^^\ 



tence of saddle node bifurcations (responsible for the cre- 
ation of the stable solutions) and also bifurcations with 
only unstable solutions. A statistical approach to esti- 
mate the size of the basin of attraction of the stable so- 
lutions was performed and we were able to observe that 
phase locking values of sin closer to (for large val- 
ues of K) present the largest basins. 



Finally we studied a system with a link coupling 
strength a varying from zero (free chain) to one (ring) to 
investigate the birth of these solutions. We have observed 
two basic processes responsible for the generation of ring 
solutions from the open chain: deformations that creates 
hysteresis for a finite range of a; spontaneous creation 
that either creates solutions coming from infinity or gen- 
erates closed circuits with four solutions (only one being 
stable) . 
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